II. Soving the diffusion problem with implicit Euler integration & finite differences


In [255]:
using PyPlot
PyPlot.version


Out[255]:
v"1.3.1"

Problem description

For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.

\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}

We therefore want to solve the following PDE:

$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$

For this second solution, we use backward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:

$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^{n+1} $$

Here, A is the finite differences matrix for the 2D Laplace operator. However, as $u^{n+1}$ is not known in advance, we have to solve this equation for it. The scheme becomes the following, which has to be solved for $u^{n+1}$:

$$ \left( I - D \cdot \Delta t \cdot A \right) u^{n+1} = u^{n} $$

Define (rectangular) problem domain and resolution for solver


In [256]:
D = 0.01 # diffusion coefficient;

In [257]:
# define domain of the problem (rectangle) and time interval
xstart = 1
xend   = 3
ystart = 2
yend   = 5
tstart = 0
tend   = 10


Out[257]:
10

In [258]:
# define number of intervals in space
Nx = 64
Ny = 32


Out[258]:
32

In [259]:
dx = (xend - xstart) / Nx
dy = (yend - ystart) / Ny


Out[259]:
0.09375

For implicit Euler, the solution should be stable independent from the time step. We can therefore choose it freely here.


In [260]:
Nt = 1000
dt = (tend - tstart) / Nt


Out[260]:
0.01

Exact Solution

As a very simple scenario, we use the following initial conditions:

$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$

This leads to the following analytical solution:

$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$

In [261]:
C0 = 5
exact(x,y,t) = C0 * sin((x-xstart)/(xend-xstart) * pi) * sin((y-ystart)/(yend-ystart) * pi) *
    exp( - (1/(xend-xstart)^2 + 1/(yend-ystart)^2) * pi^2 * D * t)


Out[261]:
exact (generic function with 1 method)

In [262]:
C_exact = Float64[exact(x,y,tend) for x=xstart:dx:xend, y=ystart:dy:yend]
pcolormesh(C_exact)
colorbar()
clim((0,C0))


Create grid and set up initial conditions


In [263]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,0) for x=xstart:dx:xend, y=ystart:dy:yend];

Define matrix for finite differences

We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).

Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:

\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}

The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.

\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}

The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.

$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$

In [264]:
# number of unknowns in x and y direction
ux = (Nx-1);
uy = (Ny-1);

In [265]:
# finite difference matrix for ∂u/∂x
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1));
FDx = blkdiag({FDx_local for i in 1:uy}...);
# finite difference matrix for ∂u/∂y
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux));

In [266]:
FD = 1/dx^2 * FDx + 1/dy^2 * FDy;

Do integration with backward Euler method


In [267]:
u = vcat(C[2:end-1,2:end-1]...);
iterationmatrix = speye(FD) - dt*D*FD


Out[267]:
1953x1953 sparse matrix with 9577 Float64 entries:
	[1   ,    1]  =  1.22756
	[2   ,    1]  =  -0.1024
	[64  ,    1]  =  -0.0113778
	[1   ,    2]  =  -0.1024
	[2   ,    2]  =  1.22756
	[3   ,    2]  =  -0.1024
	[65  ,    2]  =  -0.0113778
	[2   ,    3]  =  -0.1024
	[3   ,    3]  =  1.22756
	[4   ,    3]  =  -0.1024
	⋮
	[1888, 1951]  =  -0.0113778
	[1950, 1951]  =  -0.1024
	[1951, 1951]  =  1.22756
	[1952, 1951]  =  -0.1024
	[1889, 1952]  =  -0.0113778
	[1951, 1952]  =  -0.1024
	[1952, 1952]  =  1.22756
	[1953, 1952]  =  -0.1024
	[1890, 1953]  =  -0.0113778
	[1952, 1953]  =  -0.1024
	[1953, 1953]  =  1.22756

In [268]:
for i in 1:Nt
    u = iterationmatrix \ u
end
C[2:end-1,2:end-1] = reshape(u,ux,uy)
pcolormesh(C)
colorbar()
clim((0,C0))


Compare to exact solution


In [269]:
vecnorm(C - C_exact)


Out[269]:
0.015927455277524333

In [275]:
# collected values for varying t (Nx = Ny = 100)
error_t = [
    0.06065398100262297  # Nt = 200
    0.03290679808189858  # Nt = 400
    0.019022503842839752 # Nt = 800
    0.012077678269085844 # Nt = 1600
    0.008604595495926204 # Nt = 3200
]
Nts = [200, 400, 800, 1600, 3200]
loglog(Nts, error_t)
grid(true)



In [276]:
# collected values for varying x (Ny = 32, Nt = 1000)
error_x = [
    0.25056651685694203  # Nx = 4
    0.09271868715123792  # Nx = 8
    0.03736606527816809  # Nx = 16
    0.01958419571420384  # Nx = 32
    0.015927455277524333 # Nx = 64
]
Nxs = [4 8 16 32 64]'
loglog(Nxs, error_x)
grid(true)